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Abstract 

Background: During mammalian pre-implantation embryonic development dramatic and orchestrated changes 
occur in gene transcription. The identification of the complete changes has not been possible until the development 
of the Next Generation Sequencing Technology. 

Results: Here we report comprehensive transcriptome dynamics of single matured bovine oocytes and pre-implantation 
embryos developed in vivo. Surprisingly, more than half of the estimated 22,000 bovine genes, 1 1,488 to 12,729 involved 
in more than 100 pathways, is expressed in oocytes and early embryos. Despite the similarity in the total numbers of 
genes expressed across stages, the nature of the expressed genes is dramatically different. A total of 2,845 genes were 
differentially expressed among different stages, of which the largest change was observed between the 4- and 8-cell 
stages, demonstrating that the bovine embryonic genome is activated at this transition. Additionally, 774 genes were 
identified as only expressed/highly enriched in particular stages of development, suggesting their stage-specific 
roles in embryogenesis. Using weighted gene co-expression network analysis, we found 12 stage-specific modules 
of co-expressed genes that can be used to represent the corresponding stage of development. Furthermore, we 
identified conserved key members (or hub genes) of the bovine expressed gene networks. Their vast association 
with other embryonic genes suggests that they may have important regulatory roles in embryo development; 
yet, the majority of the hub genes are relatively unknown/under-studied in embryos. We also conducted the first 
comparison of embryonic expression profiles across three mammalian species, human, mouse and bovine, for 
which RNA-seq data are available. We found that the three species share more maternally deposited genes than 
embryonic genome activated genes. More importantly, there are more similarities in embryonic transcriptomes 
between bovine and humans than between humans and mice, demonstrating that bovine embryos are better 
models for human embryonic development. 

Conclusions: This study provides a comprehensive examination of gene activities in bovine embryos and 
identified little-known potential master regulators of pre-implantation development. 
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Background 

Mammalian pre-implantation embryonic development 
is a complex process including fertilization, cleavage 
divisions, compaction, and blastulation. During this pro- 
cess, massive degradation of oocyte-stored maternal RNA/ 
proteins and gradual activation of the embryonic genome 
take place [1]. Earlier studies employing RNA polymerase 
II inhibitor suggested that the timing of EGA is correlated 
with the speed of embryonic development. For example, 
a-amanitin halted embryo development at the 2-cell stage 
in mice [2-4] and between 4- and 8-cell stages in humans 
[5]. However, the exact timing of EGA in bovine is still 
debated. Developmental block of cultured bovine embryos 
occur between the 8- and 16-cell stages [6-8], suggesting 
EGA at this transition. Similar conclusion was reached by 
studying expression profiles of bovine in vitro embryos 
using microarray [9] or the RNA sequencing (RNA-seq) 
technology [10]. However, a microarray study utilizing 
pooled in vivo bovine embryos suggested that bovine EGA 
occurs between the 4- and 8-cell stages [11]. 

To date, the most comprehensive transcriptome profil- 
ing of bovine in vivo embryos was carried out using the 
Affymetrix GeneChip Bovine Genome Array. Although 
this microarray contains roughly 23,000 bovine transcripts, 
these represent only 12,752 genes (personal commu- 
nications with Affymetrix), approximately half of the 
mammalian genes. Previous data are also impacted by 
hybridization variations of microarray and the use of 
pooled embryos [8,11]. Although more comprehensive 
data were obtained using RNA-seq, they were restricted 
to bovine blastocysts only [12,13] or with the use of in vitro 
embryos [10]. As a result, the complete descriptions of gene 
activities during bovine in vivo embryonic development 
are still not available. Moreover, the timing of EGA should 
not be established using expression data from half of the 
genome or in vitro embryos. The RNA-seq technology 
provides unique benefits for studying gene expression 
with high resolutions and reproducibility, as well as for 
detecting novel transcripts and alternative splicing events 
[14,15]. Here we applied the Solid RNA-seq platform on 
single in vivo matured oocytes and in vivo developed 
embryos from the 2-cell to the blastocyst stages and 
obtained their comprehensive transcriptome dynamics. 
The identification of highly connected, yet relatively 
little known or completely unknown hub genes that are 
potentially master regulators of gene expression opens 
up unprecedented opportunities for further understanding 
of early development. Furthermore, we used RNA-seq 
datasets recently generated in humans and mice and 
carried out a comprehensive stage-specific comparison 
across the three mammalian species. We found that the 
three species shared more maternally deposited genes 
than embryonic genome activated genes. More significantiy, 
there are more similarities between bovine and human 



embryonic transcriptomes than those between humans 
and mice. The data obtained here will function as a 
comprehensive reference base for embryos generated 
from reproductive biotechnologies in the bovine as well 
as in the human for which the use of in vivo embryos is 
highly limited. 

Results 

Expression profiles of bovine in vivo matured oocytes 
and pre-implantation embryos 

Using linearly amplified RNA from single oocytes/embryos, 
we obtained approximately 430 million sequencing reads 
from duplicate samples of bovine in vivo pre-implantation 
embryos at 8 stages of development (Additional file 1: 
Table SI). The raw FASTQ files and normalized read 
counts per gene are available at Gene Expression Omnibus 
(GEO) (www.ncbi.nlm.nih.gov/geo) under the accession 
number GSE59186. High Pearson correlation coefficients 
were found among biological replicates of the same de- 
velopmental stage (Figure lA, Additional file 2: Table S2), 
demonstrating the reproducibility of sample preparation 
and the sequencing technology. Instead of hundreds of 
expressed genes reported earlier [8,11], the total numbers 
of detectable genes ranged from 11,488 to 12,729 from oo- 
cytes to blastocysts in our study with the use of RNA-seq 
(Table 1 and Additional file 3: Table S3). For the first time, 
it is revealed that the bovine oocytes and early embryos 
expressed roughly 50% of the total estimated 22,000 genes 
in the bovine genome. 

Pearson correlation coefficients and principal compo- 
nent analyses (PCA) of all detected genes revealed two 
distinct segmentations of bovine pre-implantation devel- 
opment: the first from the oocyte to the 4-cell stage; 
and the second from the 8-cell to the blastocyst stage 
(Figure lA, B). This segmentation demonstrated that 
EGA in bovine occurs between the 4- and 8-cell stages. 
This timing concurred with the conclusion by Kues et al. 
[11] yet contrasted with those of all other studies [6-8,10]. 
Two additional segmentations were also worth-noting: the 
first from oocyte to 2-/4-cell stages, and the second from 
8-cell/morula to the blastocyst stage (Figure IB). These 
divisions were likely results of dramatic degradation of 
maternal RNAs and early differentiation, respectively. 

Differentially expressed genes during bovine in vivo 
pre-implantation development 

Although the total numbers of genes expressed by embryos 
of different stages did not vary much, the actual genes 
expressed during early development are dramatically dif- 
ferent. A total of 2,845 unique genes were identified to be 
differentially expressed between all consecutive stages of 
development (P < 0.05). Similar to Pearson correlations on 
all detected genes, hierarchal clustering of the differen- 
tially expressed genes also partitioned pre-implantation 
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Figure 1 Correlation and hierarchical analyses of transcriptomes of bovine in vivo developed oocytes and embryos. (A). Heatmap of 
duplicate samples of the same stages of bovine embryonic development The color spectrum, ranging from red through white to blue, represents 
Pearson correlation coefficients ranging from 1 to 0.53, indicating high to low correlations. All duplicate samples are highly correlated in Pearson 
coefficients demonstrating the reproducibility of the procedures. (B). Principal component analysis (PCA) of the transcriptomes for seven stages of 
in vivo developed bovine embryos and oocytes. Embryos from the same stage are shown by symbols of the same shape. The arrows indicate the 
developmental direction of the embryos. PCI, PC2 and PC3 represent the top three dimensions of the differentially expressed genes among the 
preimplantation embryos. (C). Hierarchical clustering of differentially expressed genes in in vivo developed bovine oocytes and embryos. Two major 
clusters are shown, one consisted of the matured oocytes and embryos at the 2- and 4-cell stages. The second cluster is consisted of embryos at the 
8-, 16-cell, early morula, compact morula and blastocyst stages. The clear separation of embryos into two groups demonstrated the timing of EGA 
in cattle at the 4- to 8-cell transition. The color spectrum, ranging from red through yellow to blue, indicates normalized levels of gene expression. 
(D). The numbers of differentially expressed genes in consecutive stages of bovine in vivo pre-implantation development (P < 0.05). 



Table 1 The numbers 
matured oocytes and 



of genes detected in bovine in vivo 
each stage of in vivo embryonic 



development 


Stage 


No. of genes (FPKIVl > 0.1) 


Oocyte 


11488 


2-cell 


12678 


4-cell 


12718 


8-cell 


12729 


16-cell 


11973 


Early morula 


11670 


Compact morula 


11910 


Blastocyst 


11924 



FPKM: fragments per kilobase of exon per million fragments mapped. 



development into two distinct clusters (Figure IC), that 
from oocytes to 4-cell and from 8-cell to blastocysts, 
which confirms the timing of EGA. The majority of the 
differentially expressed genes, 2,031, were found between 
the 4- to 8-cell stages, providing another confirmation that 
bovine EGA occurs at this transition. Among these genes, 
1,086 and 945 were down- and up-regulated, respectively 
(Figure ID). The down-regulated genes are involved in 
reproduction, transcription and cell cycle regulations. 
Conversely, the up-regulated genes, representing those 
transcribed from the embryonic genome, are involved 
in translation, ATP metabolic process as well as RNA 
processing (Additional file 4: Table S4). The second 
largest change in gene expression occurred from compact 
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morula to blastocyst, with 829 genes differentially expressed 
(Figure ID), suggesting that specific genes were necessary 
during the blastulation and early differentiation processes. 
The biological processes significantly represented at this 
transition included cell proliferation, transport and early 
differentiation (Additional file 5: Table S5). This burst of 
protein production and cell division may be necessary to 
prepare the blastocyst for the upcoming coordinated 
differentiation. 

Two minor EGA events were also identified in 
addition to the major EGA between the 4- and 8-cell 
stages: one from the oocyte to the 2-cell stage and the 
other from the 16-cell to the early morula stage, with 
324 and 413 genes differentially expressed, respectively 
(Figure ID). Between the oocytes and the 2-cell embryos, 
166 of the 324 differentially expressed genes were down- 
regulated. These represent rapid degradation of the 
maternally stored transcripts. Gene ontology analysis 
indicated significant over-representation of elements 
involved in cell cycle and mitosis II (Additional file 6: 
Table S6), suggesting that the 2-cell embryos repro- 
grammed its cell cycle regulation from that of an 
arrested state to an active mode of cell division through 
changes of gene expression. The second minor EGA, 
between the 16-cell embryo and early morula, included 
413 differentially expressed genes (Figure ID). These 
genes may play important roles during the development 
of tight junctions and other processes of compaction. 
Interestingly at this transition, we found a high enrich- 
ment of genes involved in stem cell maintenance and 
development, suggesting that genes for pluripotency are 
active long before the formation of the inner cell mass 
(Additional file 7: Table S7). 

In spite of the aforementioned differences, there are 
also wide-spread similarities in the expression profiles 
between the 2- and 4-cell embryos, the 8- and 16-cell 
embryos, as well as the early and compact morulae 
(Figure ID). Specifically, only 97 genes were differentially 
expressed between the 2- and 4-cell embryos (Additional 
file 8: Table S8). Moreover, among the more than 11,000 
genes commonly expressed by both the 8- and 16-cell 
embryos, only 236 genes were differentially expressed 
(Additional file 9: Table S9). Likewise, as few as 187 
differentially expressed genes were found among the 
10,843 commonly expressed genes between the early 
and compact morulae (Additional file 10: Table SIO). 

To confirm the throughput results from RNA-seq, we 
performed quantitative real-time PGR (qRT-PCR) on 10 
genes using in vivo embryos at the 4- and 8-cell stages 
(n = 3). Among the selected genes, five {GATA6, GNB2L1, 
BAD, H2AFZ and NANOG) were up-regulated and five 
{GDP9, DNMTl, ZP2, STATS and OOER) were down- 
regulated between these two stages. The qRT-PCR results 
substantiated those from RNA-seq (Table 2). 



Table 2 Quantitative real-time RT-PCR (qRT-PCR) results of 
10 selected genes between 4- and 8-cell stage embryos 



Comparison 


Gene symbol 


Log 
(fold change) 
RNA-Seq 


Expression 


Log (fold 
change)* 
qRT-PCR 




GATA6 


8.4 


Up 


10.6 




GNB2L1 


8.2 


Up 


9.2 




BAD 


6.4 


Up 


84 




H2AFZ 


8.2 


Up 


10.3 


4- vs. 8-cell 


NANOG 


11.5 


Up 


7.8 


embryos 


GDP9 


-3.8 


Down 


-2.1 




DNMTl 


-3.2 


Down 


-2.6 




ZP2 


-3.5 


Down 


-4.4 




STAT3 


-2.8 


Down 


-2.3 




OOER 


-2.7 


Down 


-2.0 



*Fold change is expressed as the ratios of the values of the 4-cell embryos 
(n = 3) divided by those of the 8-cell embryos (n = 3}. Real time RT-PCR results 
substantiated the differential gene expression patterns from RNA-seq. 



Cluster profiles of differentially expressed genes 

Although as many as 2,845 differentially expressed genes 
were identified, the pattern of their dynamic changes can 
be categorized into as few as 30 distinct clusters , labeled as 
Clusters 1 to 30 (Additional file 11: Figure SI; Additional 
file 12: Table Sll). These clusters can be further assigned to 
four main groups of different dynamic patterns (Figure 2A). 
The first group, including Clusters 1, 7, 9, 11, 13, 14, 15, 17, 
19, 22, 24, 27 and 29, represents genes that increased their 
expression levels during pre-implantation. All clusters in 
this group, with the exception of Clusters 11, 17 and 29 
showed a dramatic increase at the 8-cell stage, indicating 
that they are transcribed from the embryonic genome. 
Genes in these clusters include developmentally important 
ones such as GATA6, H2AFZ and NANOG. Interestingly, 
genes in Clusters 11 and 29, including GATA3 and DSP, 
peaked at the blastocyst stage, suggesting their roles in 
blastocyst formation and early differentiation. The second 
dynamic expression pattern, including Clusters 2, 5, 6, 8, 
16, 18, 23, 25, 26 and 28, represents genes that underwent 
an overall trend of decrease, suggesting continued degrad- 
ation over the pre-implantation period. Of special interest 
was a sharp decrease during the 4- to 8- cell transition in 
Clusters 6, 8 and 28, including oocyte-specific genes such 
as ZP2 and WEE2, demonstrating the lack of their involve- 
ment in embryo development. The decrease of genes in this 
group may also be a pre-requisite for EGA. The third 
expression pattern, including Clusters 4, 10, 20 and 21, 
contains genes that first increased and then decreased 
their expression levels. Among these. Clusters 4 and 10 
were up at the 8-cell stages and then declined, suggesting 
that these genes, such as BAD, APOPTl and GNATl, are 
only involved in the activation of the embryonic genome. 
The last group, including Clusters 3, 12 and 30, represents 
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Figure 2 Distinct patterns and dynamics of gene expression during bovine in vivo pre-implantation development. (A) Representative 
clusters of expression dynamics during early development. Genes were clustered to be increased (a), decreased (b), increased first and then 
decreased (c), and maintained relatively constant levels of expression (d). (B). Identification of stage-specific/enriched genes by cluster analysis. 
Groups of genes were found to be only expressed in oocytes and blastocysts, enriched in oocytes to 4-cell embryos, and 8-cell to blastocysts. 
The color spectrum, ranging from red to white, represents high to low levels of gene expression. 



genes that maintained relatively constant levels of expres- 
sion throughout all stages studied. Members from this 
group such as ATPlAl, ATP5F1 and RALB are involved in 
ion exchange, energy metabolism and signal transduction, 
suggesting their necessary roles during the entire process 
of pre-implantation development. It is possible that some 
of the transcripts in this group were stable maternal RNAs 
that were never degraded while others were maintained 
by the early embryos. Nonetheless, genes in this group 
are good loading control candidates for gene expression 
quantifications. 

In addition to the dynamic changes of gene expression, 
we also identified genes that are only enriched in one 
particular stage of development. Specifically, a group of 
119 genes were enriched only in the matured oocytes 
(Figure 2B, Additional file 13: Table S12). These transcripts 
were degraded after fertilization and remained suppressed 
during embryo development. Apart from well-known/ 
studied genes such as HIFOO, we identified many less 
known/not annotated genes such as LOC78217S and 
LOC536606. These genes likely have limited roles in 
embryo development but are important in maintaining 
oocytes at the matured stage. Further investigations into 
their roles in oocytes will enhance our understanding of 
the mechanism for meiotic arrest. Another group of 234 
genes, including DNMT3A, GATA3, CD9 and APOPl, 
were only enriched at the blastocyst stage (Figure 2B, 
Additional file 13: Table S12). Groups of genes were also 
found to be enriched in a short duration of development 
such as 310 during oocyte to 4-cell stage and 111 during 
8-cell to blastocyst stage (Figure 2B, Additional file 13: 
Table S12). 



Stage-specific and cross-species gene expression 
comparisons 

In addition to analyzing changes of individual genes, we 
also examined gene-interactions by identifying modules 
of genes that were co-expressed. Gene co-expression 
suggests their involvement in a common network of 
biological processes and functions [16]. Using weighted 
gene co-expression network analysis (WGCNA) [17,18] 
we identified 17 distinct co-expression modules from 
13,127 detected genes in our RNA-seq dataset (Figure 3A). 
Twelve of these modules were stage-specific, i.e., these 
modules included genes that were overexpressed in a 
particular embryonic stage (Figure 3B, Additional file 14: 
Table S13) and can be used to represent the correspond- 
ing stage of development. Interestingly, analysis of the 
functions of genes in these modules revealed a sequential 
progression of stage-specific core gene networks. It 
migrated from cell cycle in oocytes, to regulation of 
transcription in 4-cell embryos, to translation in 8-cell 
embryos, to stem cell development, maintenance and 
differentiation in morulae, and finally to cell proliferation 
and protein transport in blastocysts (Figure 3C). Such 
coordinated changes of functional pathways are reflective 
of the little-known developmental programming. 

To explore the conservation and divergence of genes 
in the 12 stage-specific co-expression modules within and 
across species, we downloaded the raw datasets from two 
previously published microarray studies of the bovine 
[9,11] and one recently published RNA-seq study of the 
human and mouse [19]. We then identified 8,103, 9,648 
and 8,705 commonly expressed orthologs from bovine, 
humans and mice against our expression dataset. The 
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Figure 3 Co-expression network analysis of bovine pre-implantation development. (A). Hierarchical cluster tree showing modules of 
co-expressed genes identified by WGCNA, A total of 17 co-expressed modules were found during bovine pre-implantation development and were 
represented by branches and labeled by different colors to the right of the tree. The height (X-axis) indicates levels of correlation. (B), Heatmap of 
correlations (and corresponding P-values) between co-expressed modules and stage of development The color scheme, from blue through white to 
red, indicates the levels of correlation, from low to high. The stage-specific modules identified are highly correlated (I.e. over-expressed) with distinct 
developmental stages (columns). (C). Functional terms of stage-specific modules of co-expressed genes during bovine pre-implantation development 
A systematic and sequential changes in functions of co-expressed genes were observed as the embr/os progress through early development 



modulePreservation function of WGCNA was used to 
calculate the Z-statistics [20], which is a measure of the 
level and pattern of the connectivity of co-expressed 
genes. As expected, the 12 stage-specific modules were 
more significantly preserved within species than between 
species. Specifically, 5 out of the 12 bovine stage-specific 
modules were strongly (Z > 10) preserved in the two 
published bovine microarray datasets of similar samples, 
another 5 were weak to moderately preserved (2 < Z < 10; 
Figure 4). To date, there is only one published report 
[19] of cross-species comparisons using co-expression 
module data from human and mouse embryos. Here we 
conducted the first study assessing the cross-species 
preservation using data from the three available mamma- 
lian species, bovine, human and mouse. Remarkably, most 
bovine stage-specific modules were at least moderately 
preserved with those of the human but less with those of 
the mouse (Figure 4), suggesting that the human not only 



share more similarities with the bovine than with the 
mouse in genome sequences [21,22] but also in embryonic 
gene-expression patterns, and thus supporting the notion 
that bovine embryos are better models for human embry- 
onic development than their mouse counterparts. 

To further characterize the conservation and variation 
of functional modules among species, we conducted 
module analysis using WGCNA on all detected genes, 
14,766 and 13,879, respectively, from the RNA-seq data- 
sets of the human and mouse [19]. We then compared 
the gene lists within the co-expressed modules of the 
three species. Intriguingly, there are significant overlaps of 
genes in the bovine and human stage-specific modules 
(P < 10 *; Figure 5A, Additional file 15: Table S14). Of 
note, many genes in modules specific to the bovine 
oocyte, 4-cell, 8-cell and morula stages overlapped with 
those at the corresponding stages in humans. For example, 
513 genes (P < 10 '^°) overlapped between bovine and 
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Figure 4 Heatmap of module preservation of stage-specific gene 
co-expression among bovine, human and mouse oocytes and 
embryos. Commonly expressed orthologs from the present study and 
those published previously as indicated on the X-axis, including two 
from bovine microarray studies [9,1 1] and an RNA-seq study in humans 
and mice [19] were identified and included. The labels on the Y-axis 
are stage-specific modules of co-expressed genes. The color scheme, 
from white to red, indicates low to high levels of preservation. 



human oocytes alone. These genes are involved in protein 
transport and cell cycle processes. Similarly, highly sig- 
nificant overlap in module genes were observed during 
bovine and human late pre-implantation development 
(from 8-cell to blastocyst) although some bovine stage- 
specific module genes can be found in multiple stages of 
human development, and vice versa (Figure 5A). GO 
analysis of these overlapped module genes indicate 
significant over-representation of translation, RNA pro- 
cessing, generation of precursor metabolites and energy. 
To our surprise, as many as 95 genes (P < 10"^) in the 
bovine oocyte-specific module were found in the human 
4-cell specific module, giving the apparent suggestion 
that the progression of embryo development in the 
bovine may be more rapid than that of humans. This 
certainly contradicts with the established observations that 
bovine in vivo embryo development is 8 days (oocytes to 
blastocysts) [23] while in humans it is 5 days [24]. One 
possible explanation to this is the diversity of embryonic 
programming, e.g., human embryos prepare for an inva- 
sive type of implantation while bovine embryos only 
attach to the uterus [25]. Consistent with this possibility, 
GO analysis showed that these overlapped genes are 
related to signal transduction such as Ras and small 
GTPase signaling. Together, these results suggested that 
bovine and humans share many core transcriptional 
programming, while differ in stage-specificity and timing. 



In the contrary, overlaps of genes between bovine and 
mouse stage-specific modules were only observed before 
EGA and after morula formation (Figure 5B, Additional 
file 15: Table S14). For example, genes in the mouse 
oocyte-specific module overlapped significantly with those 
specific to bovine oocyte (P < 10"^"), 2-cell (P < 10"^) and 
4-cell stages (P < 10"*). Meanwhile, genes specific to the 
mouse morula module significantly overlapped with those 
in the bovine 8-cell (P < 10"^), compact morula (P < 10"^) 
and blastocyst modules (P < 10"^^). These results showed 
that mouse early and late pre-implantation genes are 
spread over a large period of the bovine development, 
consistent with the speed of embryo development in these 
two species. 

Collectively, our results show that the three mammalian 
species share more maternally deposited genes than those 
expressed after EGA. Based on overlapped genes found in 
modules prior to EGA, bovine maternal detritus (RNA 
and protein) occurs later than that in the mouse but 
slightly earlier than that in the human, despite the longer 
bovine pre-implantation development. 

Identification, visualization and validation of hub genes 

In order to identify genes that are central and highly- 
connected within the stage-specific modules, we conducted 
hub gene identification analysis. Hub genes are highly 
correlated within the stage-specific modules and are 
conceptual and concrete representatives of the corre- 
sponding modules. For each stage-specific module, we 
assigned all genes with Pearson correlation coefficients 
greater than 0.9 as its hub genes (Table 3). Furthermore, to 
explore the connections among hub genes, we examined 
the top 200 connections of the top 150 hub genes (highly 
correlated hub genes) for each stage-specific module and 
visualized them in VisANT (Table 4, Figure 6). The full lists 
of these genes can be found in Additional file 16: Table S15. 
Although there are well-studied genes such as RALB in 
oocytes and DNMT3A in blastocysts, many of these 
genes are surprisingly either under-studied in embryonic 
development or un-annotated, and are thus less known/ 
unknown in this process. For example, LOC100137763 
and LOC100849216 were highly correlated, un-annotated 
hub genes in bovine mature oocytes and blastocysts, 
respectively (Table 4, Additional file 17: Figure S2). The 
highly correlated hub genes reported here are likely key 
players for their specific stage (s) of development and 
may function as "master regulators" of gene expression and 
stage transition in early development. Further investigation 
into their identities and functions will greatly enhance our 
understanding of embryo development and our ability to 
manipulate embryos through biotechnologies. 

Cross-species analysis showed higher degrees of hub 
gene validation at the oocyte and blastocyst stages than 
at other stages. For example, 32% (211 genes) and 48% 
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Figure 5 Divergence of stage-specific gene co-expression among bovine, human and mouse oocytes and embryos. (A) Heatmap of 
gene overlap between independently constructed bovine and human modules. The X- and Y-axes show human (n = 9) and bovine stage-specific 
modules (n = 12), respectively. Each cell contains the number of intersecting genes and the corresponding P-value (-log 10) of the intersection. 
Significant gene overlaps were found in nearly all stages of development between the human and bovine. (B). Heatmap of gene overlap between 
independently constructed bovine and mouse modules. The X- and Y-axes show mouse (n = 9) and bovine stage-specific modules (n = 12), 
respectively. Each cell contains the number of intersecting genes and the corresponding P-value (-loglO) of the intersection. Significant gene 
overlaps were only observed in oocytes and morula/blastocysts between the bovine and mouse. 
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Table 3 The numbers and validation results of hub genes in bovine in vivo oocytes and embryos 



Stage-specific modules 


Total no. of 
genes/module 


Total no. of hub 
genes/module 


No. (%) hub genes (validated 
in at least one dataset)* 


No. (%) hub genes (validated 
in Kues et al. 2008 [11]) 


Oocyte_1 


2347 


650 


211 (32%) 


117 (18%) 


Oocyte_2 


815 


125 


24 (1 9%) 


6 (5%) 


2-cell 


444 


112 


12 (11%) 


11 (10%) 


4-cell 


1868 


299 


54 (1 8%) 


1 9 (6%) 


8-cell 1 


229 


52 


7 (1 3%) 


5 (10%) 


8-cell_2 


640 


172 


14 (8%) 


13 (8%) 


16-cell 


247 


27 


1 (4%) 


0(0) 


Early morula 


120 


15 


2 (1 3%) 


2 (13%) 


Compact morula 


311 


41 


7 (1 7%) 


7 (17%) 


Blastocyst_1 


1049 


274 


132 (48%) 


116 (42%) 


Blastocyst_2 


1366 


366 


59 (16%) 


44 (1 2%) 


Blastocyst_3 


1010 


118 


65 (55%) 


44 (37%) 



More than one modules were found for oocytes, 8-cell embryos and blastocysts. 

*two bovine microarray [9,1 1], one human and one mouse RNA-seq datasets [19] were used. 



(132 genes) of all hub genes in the oocyte_l and blasto- 
cyst_l modules, respectively, were validated in at least 
one dataset (Table 3, Additional file 16: Table S15). Fewer 
hub genes from the 2-cell to morula stage were successfully 
validated among species. For example, only 4% and 8% of 
hub genes were validated at the 16-cell and 8-cell stages, 
respectively. This low degree of validation reflected the 
divergence in the stage-specificity and timing of transcrip- 
tional programming. Unexpectedly, we observed relatively 
low number/percentage of validated hub genes against the 
two previously published bovine datasets [9,11] (Table 3), 
likely because of the low coverage of microarray and/or 
the relatively low resolution of microarray data. 

Table 4 Highly correlated hub genes in bovine 
stage-specific modules 



Stage-specific Hub genes 
module 



Oocyte_1 


5RPX, NAA30 


Oocyte_2 


LOCI 001 37763, PAX3, RALB, SMCIB, UNC13C VANGLl 


2-cell 


CAPRIN2, LACCl, LOC616167, NLRP9, IGLPl, POL 


4-cell 


CNTNAP2, TPM3 


8-cell 


L0C5 19952, LOC789391, LYSMD3, TEXAS 1, THAP8 


16-cell 


ARLIO, FAM84B, LOC790411, CCDC39 


Early morula 


LGALS9 STAC LOC100M0626 


Compact morula 


APOBR, GALNTLl, LRPS, PCDHW, RGS20, 




HOXAll, LOC781048 


Blastocyst_1 


DNMT3A, ATP6V0A4, FAM115C LGALSl, SLC9A3R1 


Blastocyst_2 


BCAM, BPlFAl, LOCI 008492 16, PU<NA3, 




SHR00M2 SLC16A7 


Blastocyst_3 


EEFZ RPLIOA, RPL38 



Multiple modules exist for oocytes and blastocysts. 



Pathways in stage-specific modules during bovine 
pre-implantation development 

Pathway analysis revealed essential signaling and meta- 
bolic networks in embryonic development. We found 
more than 100 pathways involved in a sequential order 
relative to bovine pre-implantation development, most of 
which were represented in oocytes, major EGA transition 
(4-cell to 8-cell) and blastocysts (Additional file 18: Table 
S16). Components of cell cycle, RNA degradation and 
progesterone-mediated oocyte maturation pathways 
were highly enriched before the 4-cell stage, while 
ribosome, spliceosome and proteasome pathways were 
highly represented after the 8-cell stage. Interestingly, 
pathways for oxidative phosphorylation, glycolysis, pyruvate 
metabolism, pentose phosphate and the citrate cycle (TCA 
cycle), which are critical not only for cell proliferation, but 
also for maintenance of pluripotency [26], were uniquely 
found in blastocysts. Additionally, many well-known path- 
ways including MAPK, insulin, ErbB, Wnt, mTOR and 
TGF-beta signaling were operative in bovine before the 
8-cell stage. It is noteworthy that the most prominent 
changes in biological networks occurred from oocyte to 
4-cell stage and blastocyst stage reflecting major functional 
transitions. 

Discussion 

The development of RNA sequencing technologies permits 
the study of gene regulation at an unprecedented level. 
Here, we provide the first comprehensive description of 
gene activities during in vivo bovine embryonic develop- 
ment. Most such studies had been conducted in the mouse 
[3,4,9,19,27]. However, mouse data have limited utility in 
human embryogenesis due to the large differences in gene 
expression and genome sequences as shown here and in 
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earlier studies. It is therefore important to establish the 
full expression profile database from an alternative species. 
To date, all expression profile studies using bovine embryos 
were either conducted on in vitro embryos and/or using 
the microarray [8-13,28]. The few studies employing the 
RNA-seq technology involved blastocyst stage only 
[12,13,28] except for a recent RNA-seq study using 
in vitro embryos of multiple stages [10]. None of these 
reports, however, can be used as the complete "gold 
standards" for bovine embryo development because in vitro 
developed embryos have wide-spread gene expression 
anomalies and the DNA microarray technology limits 
gene detection to only those printed [8,11] in addition to 
variations from hybridization. In this study, we applied the 
RNA-seq technology and revealed the transcriptomes of 
bovine in vivo pre-implantation development in a very 
high-throughput and quantitative manner [14]. For the 
first time, the bovine matured oocytes and early embryos 
were shown to transcribe more than half of all bovine 
genes [22]. 

The timing of EGA in bovine has long been accepted 
to be between the 8- and 16-cell stages [6-8]. Using 
in vivo embryos and microarray containing approximately 
half of the bovine genome, Kues et al. proposed a new 
timing: between the 4- to 8 -cell stage when the largest 
number of differentially expressed genes were found [11]. 
This result was confirmed in our study, also employing 
in vivo embryos but with all bovine expressed genes, and a 
more powerful throughput technology, the RNA-seq. 
However, two prior expression profile studies both using 
in vitro bovine embryos [8,10], maintained EGA at the 
8- to 16-cell transition. It has been shown through the 
use of RT-PCR that in vitro vs. in vivo embryos have 



step-wise differences in mRNA expression [29-31]. The 
difference in EGA timing among the aforementioned 
studies therefore further demonstrated that in vitro 
embryos are not suitable for establishing reference base 
of early development [32,33]. 

Another important contribution of this study was the 
discovery of patterns of gene expression and their correl- 
ation to milestones of embryo development. Four waves 
of transcriptional changes, between oocyte and 2-cell, 
between 4- and 8-cell, between 16-cell to early morula, 
and between compact morula to blastocyst, were each 
correlated to degradation of maternal RNA, major EGA, 
compaction and blastulation. These, together with the 
identification of transient, stage-exclusive gene expression, 
provide directions of future research in embryogenesis. 

Also for the first time, we identified a number of stage- 
specific modules in bovine pre-implantation development. 
They not only represent the corresponding stage of 
embryogenesis, but reveal an interesting progression of 
core gene networks from cell cycle (oocyte), to regulation 
of transcription (4-cell), translation (8-cell), stem cell 
development, maintenance and differentiation (morula), 
and finally to cell proliferation and protein transport 
(blastocyst). The identification of these orchestrated 
functional changes is among the first step to unveil the 
little-known embryonic programming, and is important 
in enhancing our ability to improve assisted embryo 
biotechnologies such as embryo culture conditions. For 
example, metabolic pathways unique to the bovine blas- 
tocysts, such as glycolysis, pyruvate metabolism and the 
pentose phosphate pathway, were identified. Their presence 
is compatible with the "Warburg effect" commonly found 
in cancer cells [34]. In this unique pattern of metabolism. 



Jiang et al. BMC Genomics 2014, 15:756 
http://www.biomedcentral.com/1471-2164/15/756 



Page 11 of 15 



glycolytic end products enter the pentose phosphate 
pathway instead of the TCA cycle [35]. Such variation 
from the somatic cells' metabolism of the TCA cycle 
[35] not only allows for rapid cell proliferation, but also 
maintains the pluripotency of the bovine blastocyst [26] . 
Using this feature of the blastocyst, specific medium 
that encourages the pentose pathway may be developed 
to increase the proportion of good embryos in the in vitro 
production system. 

Our cross-species analysis demonstrated that human 
embryos share more similarity to those of the bovine than 
the mouse in transcriptomes during early embryonic 
development. The expression profiles established in this 
report can therefore serve as a reference base for 
embryos from assisted technology from both cattle and 
humans. Interestingly, gene expression profiles unveiled 
unique developmental programming of embryos in the 
three species analyzed. At the superficial level differences 
in stage-specific modules appear to suggest that the 
bovine embryos progress slower than those of the mouse, 
but more rapid than those of the human. While this is 
consistent with the in vivo embryo development between 
the mouse (3.5 days) [27] and the bovine (8 days) but 
not between the bovine and the human (5 days) [24]. 
The inconsistency of gene expression at similar stages 
of development between humans and bovine suggest 
that the early embryos employ different pathways to 
prepare themselves for the upcoming different process 
of implantation. Moreover, our results showed that the 
three mammalian species share more maternally deposited 
genes than EGA-activated genes, concurring with the 
conclusion from a microarray study using the Bayesian 
clustering method [9] and again revealing species differ- 
ences in the programming of embryo development. 

The cellular and molecular mechanisms governing 
mammalian pre-implantation development are still poorly 
understood. Here we identified a number of hub genes 
that are critical connectors to other expressed components 
within each embryonic stage of development in the 
bovine. Their importance is "validated" by those that have 
been studied previously. For example, RALB is a highly 
correlated hub gene in oocytes and has key roles in both 
bovine [36] and Xenopus embryo development [37]. It 
is also implicated in tumorigenesis and cell proliferation in 
mice [38]. Similarly, DNMT3A is a hub gene in blastocysts. 
Studies in mice have demonstrated that DNMT3A is essen- 
tial for de novo methylation and embryo development 
[39,40]. DNMT3A is also likely essential in the bovine 
blastulation process. The observations that DNMT3A is 
significantly reduced in cloned bovine embryos [41] and 
that lower pregnancy/calving rates and abnormal develop- 
ment are commonplace in cloned fetuses are "functional 
validation" of the hub gene status for this important 
regulator of epigenetic modifications [42]. 



Most identified hub genes, however, haven't been studied 
or annotated, albeit their potential important roles in 
embryo development. For example, LOC100137763 and 
LOCI 008492 16 are highly correlated, yet un-annotated 
hub genes identified in oocytes and blastocysts, respectively. 
The hub genes identified here represent the unprecedented 
opportunities and insights offered by the RNA-seq tech- 
nology and bioinformatics. Collectively, our inventories 
of all hub genes provide a valuable resource for further 
studies of the molecular mechanisms of pre-implantation 
development. 

Although we had to conduct linear RNA amplification 
in order to yield sufficient materials from single oocytes/ 
embryos, the highly reproducible protocol employed here 
has been previously validated [43] and the correlation 
coefficient after amplification is higher than 0.94 [43]. 
Readers are cautioned, however, that the in vivo oocytes/ 
embryos used here were generated after superovulation 
treatment. Although superovulation can affect gene expres- 
sion of oocytes/embryos [44-47], it is frequently used in 
both research and production [48] because naturally ovu- 
lated/developed oocytes/embryos from single-ovulatory, 
large animals such as cattle are not very feasible. Nonethe- 
less, the ultimate "gold standards" for gene expression 
during bovine pre-implantation development can only be 
established using naturally ovulated/developed oocytes/ 
embryos. 

Conclusion 

This study provides comprehensive examinations of gene 
activities in in vivo bovine oocytes and pre-implantation 
embryos. Cross-species analysis revealed that bovine 
pre-implantation transcriptional profiles share more simi- 
larity to those of the human than the mice. The data 
presented here can be used to assess the impact of various 
assisted reproductive techniques in both bovine and 
human reproduction. 

Materials and methods 

Ethics statement 

Oocytes and embryos were obtained from healthy Holstein 
cows in the Institute of Animal Science, Xinjiang Academy 
of Animal Science, Urumqi, Xinjiang, P. R. China. The 
animal protocol was approved by the Animal Care and 
Use Committee of Xinjiang Academy of Animal Science 
(Research license 200815). 

Collection of in vivo matured oocytes and 
pre-implantation embryos 

Multiparous Holstein cows (n = 10) were synchronized 
and superovulated as described [49,50]. Artificial insem- 
ination using semen from one of three bulls with proven 
fertility was conducted at 12 and 24 hours post standing 
heat (Day 0). Donor animals were sacrificed at 30 hours. 
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and 2-4 days after estrus to collect in vivo developed oo- 
cytes and embryos at the 2- to 16-cell stages by oviductal 
flushing. Early morulae, compact morulae and blastocysts 
were collected by routine non-surgical uterine flushing on 
Days 5, 6 and 7. All oocytes and embryos were examined 
and staged under light microscopy. Only morphologically 
intact embryos meeting the standards of Grade 1 by the 
International Embryo Transfer Society were included in 
this study. Embryos were washed twice in D-PBS before 
frozen and stored individually in RNAlater (Ambion, Grand 
Island, NY) in liquid nitrogen. 

RNA isolation, linear amplification, library construction 
and sequencing 

Following the reproducible procedures of RNA extraction 
and linear amplification from our previous study [43], we 
isolated total RNA from individual oocytes/embryos using 
TRIzol (Invitrogen, Grand Island, NY) and co-precipitated 
the RNA with linear acrylamide (Ambion). The quality 
of the total RNA was examined with the Aglient RNA 
6000 Pico kit (Aglient Technologies, Santa Clara, CA) using 
the Aglient Bioanalyzer 2100. RNA was then amplified 
twice using the TargetAmp 2-round aminoallyl-aRNA 
amplification kit 1.0 (Epicentre, Madison, WI) according 
to the manufacturer's instructions. 500 ng of amplified 
RNA (aRNA) were used to construct the sequencing library 
following the manufacturer s instructions by SOLID™ Total 
RNA-seq Kit (Life Technologies, Grand Island, NY). After 
the sequencing library was prepared, we used an Agilent 
2100 bioanlyzer to analyze the quality of the libraries. The 
sequencing libraries were then barcoded, multiplexed, and 
sequenced on a 5500x1 Genetic Analyzer at the Center for 
Applied Genetics and Technology, University of Connecti- 
cut. We obtained 430 million sequencing reads with a read 
length of 75-bp from 16 single oocytes and embryos. The 
high correlation coefficients between samples of the same 
development stage demonstrated the reproducibility of the 
mediod (Additional file 2: Table S2). 

Mapping, assembly and gene expression analysis 

Sequencing adapters were trimmed using Cutadapt 
(https://code.google.eom/p/cutadapt/) and sequencing 
reads of low quality were pre-fUtered by FASTX-Toolkit 
before mapping (http://hannonlab.cshl.edu/fastx_toolkit/), 
using the options "fastq_quality_trimmer-Q 33-v-t20-l 
30-1". The quality of reads after filtering was examined 
using 'fastQC (http://www.bioinformatics.babraham.ac. 
uk/projects/fastqc/). Filtered reads were mapped to the 
Btau_4-6.1 assembly using Tophat [51]. Individual mapped 
reads were fed to Cufflinks [51] to construct transcriptome 
models. Any novel genes and transcripts that did not fit 
the supplied gene models (NCBI RefGene) were also 
assembled. Cuffmerge [51] was used to converge individ- 
ual transcriptome to produce a master gene model. Genes 



and transcripts mapped to uncertain chromosomes and 
contigs were eliminated. 

The merged gene model and mapping result BAM files 
from each RNA-seq library were used to quantify the 
expression levels of all genes by calculating the number 
of reads falling into each gene with the Python package, 
HTSeq [52]. A matrix of Pearson correlation coefficient 
was created using R, which was in turn used to create the 
heatmap. Principle component analysis (PCA) was analyzed 
by using R. Differentially expressed genes between two 
consecutive developmental stages were identified using 
default parameters in DESeq [52]. In each comparison, 
only genes whose sum of expression across all compared 
samples was greater than the 25th percentile were used. 
Genes were deemed differentially expressed between 
subsequent developmental stages if they showed a P-value 
of less than 0.05 (Negative Binomial Distribution). Expres- 
sion pattern clusters were generated by the K-means 
clustering algorithm using R. For gene expression patterns, 
correlations between pattern indicators and tested genes 
were calculated. P-values associated with correlations were 
also calculated and the Bonferroni correction was applied 
to adjust the P-value for multiple testing. Genes with 
adjusted P-value of less than 0.05 were considered to have 
followed the corresponding expression pattern. 

Detection of co-expressed gene modules 

The R package for weighted gene co-expression network 
analysis (WGCNA) [18] was used to detect co-expressed 
gene modules. A weighted gene co-expression network 
was first constructed, in which genes were nodes and 
connected with weighted edges. The connection weights 
between any pair of two nodes i and / were computed 
hy A{i,j) = (i + icorr(i,/)) where corr{i,j) is the cor- 
relation between the expression levels of nodes i and /' 
across all stages. The topological overlap matrix W 
which measures the topological similarity of every two 
genes (nodes) was calculated based on the connection 

weights as follows: W{i,i)= . i^^''^^^^''^,, where / 

° mm-\ki,kj\ + l-A{i,j) 




connect to both i and /, and = ^ ^^(': d)A{d,j) 

d 

with d indexing the nodes that connect to the node /. 

This topological overlap matrix has been shown to 
produce more biologically meaningful co-expressed gene 
clusters [17]. We computed the distance matrix Dhy D 
(«,/') = 1 - W{i,f) A dendrogram of clusters was obtained 
by applying a hierarchical clustering algorithm [53] on 
the matrix D. The dynamic cutting algorithm reported 
by Langfelder et al. [54] was used to cut the dendrogram 
to obtain the clusters of co-expressed genes. 
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An eigengene was calculated for each cluster as the 
principal component that explained the largest variance 
of the data in the cluster. It was a weighted sum of 
expression profiles of all genes in the cluster where the 
expression profile of a gene is a vector comprising the 
values of gene expression at the seven different stages. 
The eigengene served as the representative of the gene 
expression profiles in the cluster. Then, clusters whose 
eigengenes were interrelated with correlation of more 
than 0.7 were merged. The final clusters of genes were 
referred to as gene co-expression modules. 

Stage-specific module identification 

To detect modules whose eigengene showed high expres- 
sion levels at a specific stage but low in others, we used 
a unit vector to indicate each stage. In other words, the 
entry of this unit vector for the corresponding stage was 
one, and zero for the others. We then computed the 
correlations between each stage-indicator vector and 
the eigengene of each module, which also yielded a P- 
value associated with each correlation. Smaller P-values 
corresponded to more significant correlations. If a mod- 
ule received P < 0.05 for the correlation at a particular 
stage, we labeled the module as specific to that stage. 

Module preservability/reproducibility 

We downloaded and mapped genes from two microarray 
datasets of the bovine [9,11], and two RNA-seq datasets 
of the human and mouse [19], as well as our own to the 
orthologous gene database (http://www.ncbi.nlm.nih.gov/ 
homologene/). After identifying the commonly expressed 
orthologs, the preservability of a module was measured 
by the Z statistics [20], which characterizes the density 
and connectivity of genes within a module to those in 
the validation dataset. The function, module Preservation, 
in the WGCNA package was used to calculate the Z 
statistics. The categories of preservation were defined as 
strong if Z > 10, weak to moderate if 2 < Z < 10, and no 
evidence of preservation if Z < 2, as suggested by an 
early simulation study [20]. 

Cross-species module overlapping analysis 

To study if the development of functional modules 
conserves across species, we compared the gene co- 
expression modules of the bovine, mouse and human. 
The same module detection analysis was performed on 
the human and mouse datasets by Xue et al. [19]. The 
number of overlapping genes in any two modules each 
from a different species was counted. Fisher's exact test 
was conducted to show whether or not the degree of 
overlapping was simply due to a random chance, which 
yielded a P-value reflecting the statistical significance of 
the overlap. 



Module hub gene identification and validation 

The membership of a gene in a module was measured 
by the correlation between that gene and the eigengene 
of the module. Genes in a module that are highly corre- 
lated with the module eigengene are defined as hub 
genes for the module. We used all genes with correlation 
to their module eigengene of greater than 0.9 as the hub 
genes. To explore the connections among hub genes, we 
examined the top 200 connections of the top 150 hub- 
genes for each stage specific module and visualized them 
in VisANT [55]. To validate the hub genes, we used the 
raw datasets from two previously published microarray 
studies in the bovine and one RNA-seq study in the human 
and mouse. These data were subjected to WGCNA and 
stage-specific modules and lists of hub genes (IdVIE > 0.9) 
were generated for each dataset. We then determined the 
overlap of hub genes from each stage-specific module of 
the same developmental stage in different datasets. 

Gene ontology analysis 

Functional annotation enrichment analysis for Gene Ontol- 
ogy (GO) was conducted by topGo package in R [56]. Data- 
base for Annotation, Visualization and Integrated Discovery 
Bioinformatics Resource [57] was used for pathway ana- 
lyses. We summarized all similar sub-GO terms and path- 
ways into an overarching term, and P-values are shown for 
the representative terms. 

Validation of RNA-seq data 

Quantitative real-time PGR (qRT-PCR) was performed 
to validate differential expression of 10 selected genes 
using embryos at the 4- and 8-cell stages (n = 3). Among 
these, five genes {GATA6, GNB2L1, BAD, H2AFZ and 
NANOG) were up-regulated and five {GDP9, DNMTl, 
ZP2, STATS and OOER) were down-regulated between 
these two stages. Amplified RNA from individual embryos 
was reverse transcribed to cDNA by Superscript III Reverse 
Transcriptase (Invitrogen) and amplified with specific 
primers (Additional file 19: Table S17). The qRT-PCR was 
performed using SYBR Green PGR Master Mix (ABI) and 
the ABI 7500 Fast instrument. Data were analyzed using 
the 7500 software version 2.0.2 provided with the instru- 
ment. All values were normalized to the internal control, 
GAPDH. The oocytes and embryos from 2-cell to blastocyst 
stages were pooled and used as the calibrator sample. 
The relative gene expression values were calculated 
using the 2 '^^'^' method. The mean for each stage was 
determined and compared for an overall fold change. 

Additional fifes 



Additional file 1: Table SI. Summary of sequence read alignments to 
the reference genome. 
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Additional file 2: Table S2. Pearson correlation coefficients of duplicate 
bovine oocytes and embr/os of the same stage. 

Additional file 3: Table S3. Normalized read counts (expression) of 
genes in bovine oocytes and embr/os. 

Additional file 4: Table S4. Differentially expressed genes between 
the 4- and 8-cell embryos. Spreadsheet 1: all differentially expressed 
genes; Spreadsheet 2: genes down-regulated; Spreadsheet 3: genes 
up-regulated; Spreadsheet 4: GO analysis output of down-regulated 
genes; Spreadsheet 5: GO analysis output of up-regulated genes. 

Additional file 5: Table S5. Differentially expressed genes between 
compact morulae and blastocysts. Spreadsheet 1 : all differentially 
expressed genes; Spreadsheet 2: genes down-regulated; Spreadsheet 
3: genes up-regulated; Spreadsheet 4: GO analysis output of down- 
regulated genes; Spreadsheet 5: GO analysis output of up-regulated 
genes. 

Additional file 6: Table S6. Differentially expressed genes between 
oocytes and the 2-cell embryos. Spreadsheet 1: all differentially expressed 
genes; Spreadsheet 2: genes down-regulated; Spreadsheet 3: genes 
up-regulated; Spreadsheet 4: GO analysis output of down-regulated 
genes; Spreadsheet 5: GO analysis output of up-regulated genes. 

Additional file 7: Table S7. Differentially expressed genes between 
the 16-cell embryos and early morulae. Spreadsheet 1: all differentially 
expressed genes; Spreadsheet 2: genes down-regulated; Spreadsheet 
3: genes up-regulated; Spreadsheet 4: GO analysis output of 
down-regulated genes; Spreadsheet 5: GO analysis output of 
up-regulated genes. 

Additional file 8: Table S8. Differentially expressed genes between the 
2- and 4-cell embryos. Spreadsheet 1: all differentially expressed genes; 
Spreadsheet 2: genes down-regulated; Spreadsheet 3: genes up-regulated; 
Spreadsheet 4: GO analysis output of down-regulated genes; Spreadsheet 5: 
GO analysis output of up-regulated genes. 

Additional file 9: Table S9. Differentially expressed genes between 
the 8- and 1 6-cell embryos. Spreadsheet 1 : all differentially expressed 
genes; Spreadsheet 2: genes down-regulated; Spreadsheet 3: genes 
up-regulated; Spreadsheet 4: GO analysis output of down-regulated 
genes; Spreadsheet 5: GO analysis output of up-regulated genes. 

Additional file 10: Table S10. Differentially expressed genes between 
the early and compact morulae. Spreadsheet 1: all differentially expressed 
genes; Spreadsheet 2: genes down-regulated; Spreadsheet 3: genes 
up-regulated; Spreadsheet 4: GO analysis output of down-regulated 
genes; Spreadsheet 5: GO analysis output of up-regulated genes. 

Additional file 11: Figure SI. All clusters of expression dynamics 
during early bovine in vivo embryo development 

Additional file 12: Table S11. Distinct Clusters of gene expression 
patterns in bovine oocytes and embryos. 

Additional file 13: Table SI 2. Stage-specific/enriched genes in bovine 
oocytes and pre-implantation embryos. 

Additional file 14: Table S13. Co-expressed genes in stage-specific 
modules of bovine oocytes and embryos. 

Additional file 15: Table S14. Gene expression overlap between 
species. Spreadsheet 1 : overlapped genes between bovine and human 
embryos; Spreadsheet 2: overlapped genes between bovine and 
mouse embryos. 

Additional file 16: Table S15. All hub genes identified in bovine 
oocytes and pre-implantation embryos (genes in blue are validated in at 
least one dataset of the human, mouse and bovine). 

Additional file 17: Figure S2. Visualization of representative highly 
correlated hub genes in bovine oocytes and embr/os. 

Additional file 18: Table S16. Functional pathways in stage-specific 
modules in bovine oocytes and embryos. 

Additional file 19: Table S17. Primers for real time qRT-PCR. 
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